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Gravitational particle production naturally occurs during the transition from the inflationary phase to the non- 
inflationary phase. If the particles are stable and very weakly interacting, they are natural nonthermal dark 
matter candidates. We show that such nonthermal dark matter particles can produce local non-Gaussianities 
large enough to be observed by ongoing and near future experiments without being in conflict with the existing 
isocurvature bounds. Of particular interest is the fact that these particles can be observable through local non- 
Gaussianities even when they form a very small fraction of the total dark matter content. 



I. INTRODUCTION 



Standard slow-roll inflationary models with a single dynamical field degree of freedom (e.g. see the review article 
lH]) cannot generate large local non-Gaussianities (NGs) |l2K5l, which have been widely discussed and speculated upon 
in the context of the cosmic microwave background (CMB) data (e.g. 1^-12]) and large scale structure data (e.g. 1 13 - 
|22J ). Many multifield mechanisms have been proposed to generate observably large local NGs (e.g. |23-,46.|). Most 
of these models utilize coherent condensate field degrees of freedom instead of incoherent many-particle states. 

In this paper, we explore the possibility that nonthermal dark matter (DM) particles gravitationally produced during 
the phase transition out of the quasi-de-Sitter phase of inflation |47 48 1 generate observably large NGs. ' These dark 
matter particles can be viewed as the remnants of de Sitter (dS) temperature driven radiation during inflation, and no 
non-standard ingredients are needed for the inflationary scenario for the purposes of this work. The only nontrivial 
model requirement is that the dark matter either be very heavy and/or superweakly interacting. 

This class of scenarios effectively possesses only three important independent dimensionful parameters: the Hubble 
expansion rate during inflation, the reheating temperature, and the dark matter mass. Hence, the physics is dominantly 
controlled by only two of these parameters since the third converts the other two into dimensionless numbers. We 
choose these to be the dark matter mass nix and the reheating temperature Trh. The existing cosmological data 
constraining the isocurvature perturbation amplitude and the dark matter abundance place bounds on the allowed 
parametric range for these parameters. We find that to generate large observable local non-Gaussianities characterized 
by an effective f^i parameter of around 30, there is an upper bound of mx ^ 4-He, where He is the expansion rate 
at the end of inflation. We also find that f^i will be suppressed if Trh ?i 10'' GeV if the dark matter is absolutely 
stable.^ Somewhat surprisingly, even when the X particles make up a small fraction of the total dark matter content 
while thermal relics make up the rest of the dark matter, observably large non-Gaussianities may be imprinted. 

The isocurvature perturbations in this class of scenarios have been studied previously |49 |. We note that this was 
also briefly considered in |23|, which anived at a pessimistic conclusion. However, that paper did not consider the 
model as carefully as |49|, which reached a more realistic conclusion regarding the viability of such scenarios. The 
purpose of this paper is to point out that within this framework, large local non-Gaussianities can be generated with a 
single (9(10^') tuning of the dark matter mass. 

The order of presentation is as follows. In Sec. II, we discuss the class of dark matter and inflation models for which 



the current non-Gaussianity computation is relevant. In Sec. Ill we present a computation of the two-point function, 
including the cross-correlation function between the isocurvature and the curvature components. The computation of 
the bispectrum and a presentation of detailed arguments as to how the local non-Gaussianity can be large is given in 
Sec. IV In Sec.JV] we check our general analytic arguments by computing in detail numerically the observables in the 
context of m^^^l inflationary model. We then close the main body of the paper with a summary and conclusions in 
Sec. VI In Appendix [a] we present an analytic approximation to the mode function during inflation that accounts for 
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' These particles are sometimes called gravitationally produced WIMPZILLAs. 

^ This mass scale has part of its origins from the maximum dark matter abundance today. 
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the small deviation from the pure dS phase. In AppendixjB] we present the details of the computation of the curvature- 
isocurvature correlators. Finally, in Appendix |Cj we justify how the effective classical background variable about 
which the classical perturbations are defined is given by the expectation value of the quantum operator. Throughout 
this work, our metric convention is (+, — , — , — ). 

II. CLASS OF DARK MATTER AND INFLATIONARY MODELS 

We begin by defining the class of dark matter and inflationary models considered in this work. One requirement is 
that the dark matter field X be sufficiently long lived to be a viable dark matter candidate. Since we are considering an 
isocurvature dark matter scenario, another requkement is that X is very weakly interacting with thermalized Standard 
Model (SM) particles that are assumed to arise from the inflaton decay chain. Although it is straightforward to include 
dark matter interactions that allow transformations to different particles, we will assume that the self-annihilation (or 
coannihilation) interactions of the dark matter are too weak to change the dark matter number density appreciably 
after it is produced at the end of inflation. Next, we note that even if X is sufficiently weakly interacting as not 
to thermalize, it may be a byproduct of a slow-roll inflaton decay with a strength larger than gravitational strength. 
In such situations, the isocurvature nature of the X particles produced gravitationally will be made impure by the 
inflaton decay contribution. To keep the analysis simple for the purposes of this paper, we will assume that the decay 
contribution is negligible. Finally, we will assume thatX is a boson minimally coupled to gravity. We will explore the 
complexities that arise when some of these requirements are relaxed in future work. 

The simplest model that satisfies the above criteria has two scalar fields minimally coupled to gravity as follows: 

5 = J d'^x^^ [{d^f-2Ui^) + idxf-mlx^] +Ssm[8^v,Wsm}]+Srh[8^v,^,{\I/sm}1 (1) 

where Ssm is the SM sector, Smi is responsible for reheating, and is the inflaton realizing a slow-roll inflationary 
scenario with t/(0) = m^(j)^/2. We note that since particle production during the dS to non-dS phase leads to a dS 
horizon temperature population of SM one-particle states, one might naively expect a minimum reheating temperature 
of Trh ^ He/ {2k), where is the Hubble expansion rate the end of inflation. However, this is incorrect since during 
the coherent oscillations period between the end of the inflation and the radiation dominated epoch, the SM radiation 
dilutes with respect to the inflaton energy density. 

Since we will carry out numerical analysis of the mode equations, this simple model is useful. Furthermore, it is 
clear that the results generalize to a large class of models where the interactions are very small. We note also that the 
requirement of the interactions being weak enough to avoid thermalization does not require a particularly stringent 
limit on the interaction couplings. For example, for typical 0(1) coupling strengths for self-annihilations, it is well 
known that large values of mx will naturally lead to the nonthermal DM behavior desired in this paper if 

/200TeVy/7^X (2) 

where T^ax }t 7rh is the maximum effective temperature reached during reheating fSU]. Of course, one possible model- 
building obstacle with large masses is that in situations with accidental global symmetries protecting the stability of 
the particles, higher-dimension operators must be suppressed to avoid early decay. Nonetheless, many viable beyond 
SM (BSM) models that contain superheavy DM candidates have been proposed ||5TH63]| . 

Given that the X particles have negligible non-gravitational interactions and minimal couplings to gravity, they can 
only be produced gravitationally or through initial conditions. We consider the dynamics of an inflationary patch 
whose Bunch-Davies vacuum |64| satisfies 

limai|BD,0) =0, (3) 

where is the annihilation operator associated with the curvature perturbations which is approximately dominated 
by the inflaton (see Appendix [B] for more details). This vacuum is also assumed to satisfy the "no-particle" condition 
of the adiabatic vacuum during inflation ll47l I65l - l67l : 

dk\BD, 0)=0, (4) 

where d^ is the annihilation operator associated with the X field. The stress-energy tensor is renormalized such that 

(fiZ),0|f^™V£',0)=0, (5) 
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which in practice is accompHshed by the normal ordering of the creation-annihilation operators in the adiabatic vacuum 
basis. This means that the classical initial condition dependent DM density vanishes. 

Nonetheless, because the transition from quasi-dS phase of inflation to the non-dS phase after inflation represents a 
non-adiabatic transition, it is well known that non-negligible particle production occurs through Bogoliubov mixing of 
the creation-annihilation operators, giving rise to a significant DM abundance today [47 J . The physics of this particle 
production mechanism is similar to that of Hawking radiation. In the small mass limit with minimal gravitational 
couplings, we find numerically that the X energy density at the end of inflation can be approximated as Px{te) ~ 
0.01//g , which leads to the rehc abundance of X particles today to be 

axh^ « 10-' ( ,f V ( p^" ^ . (6) 

(In the actual parametric region of interest for the numerical calculations in Sec. [V] we cannot neglect the mass and 
will have a smaller relic density.) What is interesting about this scenario is that although the classical picture of particle 
production occurs at the end of inflation, the correlations that are relevant at the CMB scale are set long before the 
bulk of the particle production occurs. This is intuitively self-consistent from the Heisenberg time-energy uncertainty 
considerations. Although Eq. (|6]l yields the simplest possible scenario, we later generalize the situation to the case of 
mixed dark matter contributions in which the total cold dark matter (CDM) abundance is given by 

i^CDM = iitherm + , (7) 



where ntheim are thermal relics that have only adiabatic perturbations and £2^ are relics that have dominantly isocur- 
vature perturbations. In this case, we define 

(Ox = :^, (8) 

and will scale some of our computations to generalize our results to a wider class of scenarios. 

We now comment further on the inflationary model relevant for the above DM scenario. The main features of the 
inflationary model that are numerically important for the isocurvature and non-Gaussianity analyses are (the Hubble 
expansion rate when the modes of interest leave the horizon). He, and Trh (the reheating temperature). As we will see, 
the primary role of H^, is to determine the spectral index of the isocurvature spectrum, < //* controls the particle 
production, and Trh partially controls the map between the comoving wave vector and the physical momentum. 
We assume that there are curvature perturbations from the inflaton sector with the right magnitude to approximately 
explain the CMB spectrum. As we will see, the current observational limits require that the isocurvature contribution 
is subdominant. 

Finally, it has been noted ll23l that this class of models suffer from the boundary condition of {X) = being an 
unnatural expansion point of the fluctuations of X for nix ^ H. Although it certainly is true that in this Umit the 
H dependent radiative corrections lift the flatness of the potential, there are no radiative tadpoles that are generated. 
Furthermore, although it is true that once the non-decaying mode decoheres as the wavelength is stretched outside 
the horizon, acting like a classical background with (X) ^ over the patch of the size of that wavelength, there is 
no obvious worry that we are choosing a patch that has {X) = as long as the inflation did not last many orders of 
magnitude in efolds longer than what is needed to explain the CMB data. Indeed, because of the slight blue tilt, there 
is no infrared divergence in this class of models. As we will see, the blue tilt of significant isocurvature amplitudes is 
still compatible with observations. 

For completeness, let us also explicitly state the cutoffs of our theory. Since the spectral index of the correlator is 

related to 3 - 3 ^J\-Am\/9H'^ and our scenario has > 0, the correlator has a blue spectrum (as we will see more 
explicitly in Sec.[lll|. Thus, we do not have any IR divergences for the isocurvature correlations. The UV cutoff of our 
theory is set to be the horizon scale at the end of inflation: k\j\ — HeOe- 

Finally, we note that we use the scalar metric perturbation parameterization 

ds^ = (1 +E)dt^ - ladjFdtdx' - [5ij +A5ij + didjB] dx'dxK (9) 

We make the usual choice for the gauge invariant variable that describes the inflaton dynamical degree of freedom: 

C^t-//?^. (10) 
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In terms of this variable, the neariy scale invariant slow-roll inflaton power spectrum is given as 



where to leading order in the slow -roll parameter 



e = 



2 [ui^k) 



we have 



(11) 



(12) 



(13) 



In the above, denotes the field value when the mode k leaves the horizon. We will give more details about the X 
fluid variable when discussing the two-point function in the next section. 

In summary, the class of dark matter models that is relevant for this paper corresponds to cases with gravitationally 
produced bosonic dark matter that never fully thermalizes with the reheating radiation produced from the inflaton 
decay. The slow-roll inflationary model produces the dominant curvature perturbation spectrum and couples to the 
dark matter sector only gravitationally. 



III. TWO-POINT FUNCTION 



Although the isocurvature perturbations previously have been computed for this class of models B9I . here we redo 
the analysis with more careful attention to cross-correlations between the curvature perturbations and the isocurvature 
perturbations because the observational constraints have become increasingly severe. 

To begin, consider the energy momentum tensor of X: 



41^ ^d^XdyX-g^y 



daXd"X -V{X) 



(14) 



where V{X) = m\X^ /2. Comparing this to 



we see that if we define 



(perfect fluid) / 



(P) 
'X 



ip) 



(15) 



d^'X 



we can satisfy the equality 



if 



and 



g"PdaXdpX 



rp{X) (perfect fluid) 



p|"' = M^M^r^p;*^'""''') = ^d"xdaX+v{x) 



pl,'^ = ^d''xdaX-v{x). 



(16) 



(17) 



(18) 



(19) 



We note that to identify Eq. ( 16 1 with the fluid velocity, d^^X has to be timelike. This is consistent with the fact that 
any wave packet made of on^iEell 1 -particle states can be decomposed in terms of mode functions characterized by 

timelike 4-momenta. Unlike the coordinate dependent Tqq \ p'^'^ is a diffeomorphism scalar. We also note that even 
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though Eq. ( 18 1 looks like it has the wrong sign between the (doX)^ and | VXp, the sign is correct and p^''^ is positive 
definite whenever the fluid interpretation is valid (whenever df^X is timelike). 

We now quantize X by promoting it to an operator X. As explained in Appendix [c] this allows us to identify 

P^'^^CpI''^), (20) 

where the normal ordering is with respect to the operators defining the X vacuum state during the quasi-dS era. After 
Bogoliubov transforming : p^f^ : to operators of 1-particle states relevant for non-dS spacetime, (: j3^''' :) will develop 

a nonzero value at that later time. We note that (: pj/''' :) is homogeneous as long as the vacuum state is spatially 
translation invariant. (Here the inflaton/scalar perturbations are treated as operators, which means that as long as the 
vacuum governing these are spatially translal 
Next, we consider the semi-classical variable 



vacuum governing these are spatially translation invariant, (: p^''^ :) will be spatially translation invariant as well.) 



We can then define the usual fluid variable 



5p^/^^p^/^~p^,'\t). (21) 

Cx^f-//-^, (22) 

where the spatial scalar metric perturbation = a^(f ) [A5,y + didjB] . Under the diffeomorphism f ^ f + 1^", we have 

A A + 2H^° (23) 
5pW ^ 5p^/^+^'^^/,''\ (24) 

which makes gauge invariant, as expected. Since the gauge invariant isocurvature variable that describes the 
difference from an inflaton descendant fluid and the dark matter is conventionally defined as = 3 (i^x ^ C) : we find 
during the time period after po{t) > stops oscillating and behaves as a^^ that 

5s,^3(^+^-5x~i;), (25) 



where 



Note that we did not convert ^ in terms of the inflaton field since the inflaton eventually decays while ^ remains a 
conserved quantity on long wavelengths even after the inflaton decay if the ^ associated with the decay product fluid 
is adiabatically matched to the ^ composed of the inflaton. 

In the comoving gauge defined by the coordinate system in which the inflaton fluctuations vanish (i.e., 5(p = 0), 
C=aW/2. Hence, 

where the (c) superscript refers to the comoving gauge quantity. Therefore, correlators involving Sj^'^ and A^"^' are 

of primary physical interest. To compute them, we quantize 5p^'''''^' by promoting p^'''^'^' —t' p^f^^'^^ (through the 

quantization of X) and promoting PQ'\t) to an identity operator (since this was already defined semi-classically as a 
matrix element according to Eq. ( |C23[ )) as follows: 

5pW«=:pi^)('):-ip<'')(f). (28) 

This can be used to compute {5pl^''^^'\t,r)8p^''^^'\t,0)) and then can be divided by Po(f) after this quantity settles 
down to give an expression for {5![\t,r)5j^\t,Q)). 



k 




k 



ki — k 



Figure 1: A diagrammatic representation of the 2-point function of a composite operator Sx- The open circle corresponds to 
external momentum flowing out, and the shaded circle to external momentum flowing in. The slashes on the legs indicate that the 
X propagators are on- shell. 

Diagramatically, the 2-point function is as shown in Fig. [T] Hence, the 2-point power spectrum is 



2n^ 

{In) 



(29) 



in which 



2p, 



(30) 



5^'' is approximated as m\X^ /2Pq'\ and Xj^ is the solution to the mode equation in Appendix |a| For nix/H < 3/2 
and k/ (aH) ^ 1, we can express Px{k) approximately as 



Px{k)=Ax 



^0 



(31) 



where 



^(^)-3-3Jl-g>0. 



(32) 



In 114911 . H^, is allowed a wave vector dependence; however, this effect is a subdominant correction to the already 
small (mx/H*)^. It is important to note that for the parameter range of interest 7a: ^ 1 such that k^Px is nearly scale 
invariant, the amplitude Ax is given by the approximate formula 



Ax 



2rx( 



i/3 



- exp 



3ilrx{H,)-YxiHe) + 3ln 



(33) 



where £ 



Mi 



is the usual inflaton slow-roll parameter. The complicated exponential factor arises from 



considering the time evolution of the mode function Xj^ to a time beyond the time when k/a < e/H^,. We see that with 
a tuning of the mass parameter to 0(0. 1 ) precision. Ax can be a small number despite the exponential factor of e <C 1 . 
For the dark matter abundance to be compatible with cosmological observations, it is important that //* > He while 
mx/He ^ 0(1). Using these approximations, we find 



4 
Jx 



Px{k) 



(34) 



where we note that if k^Px{k) has a flat spectrum, will also inherit a flat spectrum. On the other hand, we will see 
in the next section that the key to obtaining large non-Gaussianities is that (k^Px(k)Y is much smafler than k^Px{k). 
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Figure 2: A diagrammatic representation of the cross-correlator {^8s), whiere tlie leading interaction vertex is proportional to mj^. 
The wavy propagator corresponds to the on-shell {^Q correlator. The open circle corresponds to external momentum flowing out, 
and the shaded circle to external momentum flowing in. The slashes on the legs indicate that the X propagators are on-shell. 



Thus far, we have focused on only the X isocurvature perturbations. In the mixed scenario described near Eq. (|7|, 
we can rescale Eq. ([34| to obtain the total CDM isocurvature perturbations as 



Ai^{k) = (o],Al (k), 



(35) 



since the rest of the dark matter contribution has no isocurvature perturbations. 

As we will see in the next section, the bispectrum is maximized in the parameter region in which A? {k) is large. 

The observational bound on A|^ is suppressed unless 1171 1691 



(C55>« V(CC)(5s5s>. (36) 
The left hand side (C^s) can be computed using the in-in formalism 0|7Ol using the trilinear interaction Hamiltonian 



H,it) 3 



3 3 

d xa 



Dimensional analysis corresponding to the dominant diagram shown in Fig. |2] tells us that 



(37) 



(38) 



Hence, this might lead to the naive conclusion that the power spectrum of { t, 5s) is approximately of the order of 

, where tg is the time at the end of the inflation. However, as shown in Appendix 

dS wave function scaling factor which is significant that is neglected in the naive estimate. Including this factor, we 
find the following cross-correlation bound: 



there is a nontrivial 



(39) 



This satisfies Eq. (36i, making the isocurvature uncorrected with curvature perturbations as long as o: > A| ^ 10 ^. 
Therefore, since the isocurvature is of the uncorrected type, the cuiTent observational bound on the adiabaticity in 
terms of the parameter a (see \T\ and references therein) is 



a < 0.067. 



(40) 



We need not woiTy about the lower bound of a > 10 ^ because the isocurvature perturbations are observationally 
irrelevant when a is that small. In the next section, we will see that this implies the maximum fyi ~ 0(50). 
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Pi, X 




Figure 3: A diagrammatic representation of the 3-point function of a composite operator Sx- The slashes indicate that the X 
propagators are on-shell. The spatial variables {x,y,z} indicate the insertion points of the external momenta p,. 



IV. BISPECTRUM 

In this section, we now compute the bispectrum Bg^{pi,p2,P3) defined by 

where we recall from the previous section that in the comoving gauge can be identified with the gauge invariant 
quantity ds^- With the diagrammatic representation given by Fig.[3j we find 

Bs,ipuP2-iPi+P2))-8(ak I ^Px{\k\)Pxi\pi+k\)Px{\P2-k\), (42) 

which gives the analytic estimate of the primordial bispectrum for different triangle configurations fixed by pi and p2- 
As large non-Gaussianities are difficult to obtain in slow-roll inflation in the squeezed triangle limit, we will focus 
on that limit in this work. Using the one-pole approximation, we can estimate the isocurvature bispectrum (written in 
a symmetrized form in the wave vectors) in the limit that one of the |/?, | is much smaller than the others as follows: 

B5siPuP2,P3) « ^^^^^Px{Pmin)(4[Px{pi)Px{p2)+Px{p2)Px{p3)+Px{p3)Px{pi)] (43) 

where py^m — min{|/5, |}. Unlike some curvaton cases, Px(Pmm) does not have an IR divergence because of the blueness 
oi Px ( Yx {H* ) > 0). Since the energy density of X is quadratic in X, the density correlator scales as and not just Px, 
which means that the coefficient P^^^J^xiPmm) is not quite the power spectrum. Nonetheless, because of the blueness 
of Px, P^^iiyPx{pmm) can be strongly suppressed if JxiH*) is large, and hence mx has to be smaller than //* for a 
non-negligible bispectrum. On the other hand, if mx is too small compared to we saw in the previous section that 
the isocurvature perturbations are larger than what is allowed by cuiTent data. Hence, there is a window for which the 
non-Gaussianities can be large and the isocurvature perturbations are consistent with the existing data. 

To see why the bispectrum composed of quadratic fields is larger than the bispectrum composed of linear fields 
(such as for ordinary inflatons), consider 

f =8x5 Bs^ipuP2,P3) 

6Pi;ipi)PdP2)+PdP2)Pdp^)+Pdp^)Pdpi)' 

where the factor of 8 comes from the fact that the ^ contribution to the temperature fluctuation AT /T compared to 5s 
contribution on long wavelengths is different by a factor of 2, i.e. 

^ = -k--5s. (45) 
T 5^ 5 
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In the squeezed triangle limit, the effective fffL is analytically estimated to be 

,3 



^160»|//^ .\ Px{pi )Px jpi) + Px {P2)PX {P3 ) + Px {pi )Px jpi ) 

^''^ 3 Yx \2n^ ""^^'""V Pi;{P^)PdP2)+PdP2)PdPi)+PdPi)PdP^)' 



(46) 



The large numerical factor of 160/3 can be traced to the factor of 8^ that arises from the product of the contraction 
permutations and the relative weight of the ^ and 5^ contributions to AT /T. As P^^^PxiPmin) is suppressed partly 
from the smallness of the Px amplitude as well as the blue tilt (causing p^^^ to be a suppressor), a large ratio of Px/Pi^ 
is required for an unsuppressed fyi- As we will now argue, a Px /Pi^ will arise from the fact that 5^ is quadratic in X. 



To understand how to obtain a large Px/Pi^, we first note that since 5^ is quadratic functional of X (see Eq. (|34|i). 



we have 



If we define lEl 



we find for a ^ 1 that 



Px 



a 

l-a 



rx^is 



1 

A- 



1 

COx' 




(47) 



(48) 



(49) 



Hence, if cOx = 1, as long as a » 10 ^/yx , there is a large ratio of Px/Pt; because A^ ^ ^% ^x- 

Combining Eq. ( [46] l and Eq. (|49]), we see 

(50) 



fNL - lO^a^/^y^. 
Hence, if we can achieve a w 0.07 Q and ^ 0. 1 in our model, we can achieve 

fNL- 0(50), 



(51) 



which is clearly observable according to the forecasts of Planck and large scale structure experiments (see e.g. 17211731 ). 
Remarkably, this is independent of (Ox- Hence, even when the dark matter composition of the X particles is small, 
non-Gaussianities associated with X may be observable. We note that hidden in this analytic estimate is the implicit 
assumption that a can remain fixed as cox 0. However, Eqs. ( [34| and ( [35| ) show perturbation theory would break 
down if a needs to be fixed as cOx ^0 limit is literally taken. For example, to have large local non-Gaussianities, the 
CDM isocurvature should be 5s — COxSx ^ a^^^i^ ~ 10^^. Hence in the case cOx ^ 10^^, the perturbativity bound 
of 5x < I is violated. In the next section, we will compute the relevant quantities more precisely in the context of a 
simple t/(0) = m50^/2 inflationary model. 



V. NUMERICAL RESULTS 



As shown in Sec. IV a local non-Gaussianity value of fy^ ^ 50 is achievable as long as the parameters of the model 
result in o: w 0.7 and Yx ^0.1. The identification of these parameters requires the computation of \Xii\^ associated 
with the mode function. Although an analytic approximation exists in Appendix |A] it is still difficult to identify the 
parametric dependence because of the fact that the massive field X, unlike the variable ^, evolves during inflation even 
when the mode wavelength is superhorizon in magnitude. Its evolution depends on the variable yx, whose slow time 
dependence is difficult to account for analytically. Hence, to check the phenomenological viability of this isocurvature 
model, we now compute the necessary mode functions numerically within a = model. 

After inflation, the energy density px (as investigated analytically and numerically in 14811741 ) is estimated to be 



Px 



(52) 
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To obtain the k that appears in the mode function Xk, we use the pivot scale ko — 0.002 Mpc and the standard 
reheating relationships i731 

(53) 
(54) 

-2/3 / \ 1/3 



k 


ak, 

= — Ko 




flo 




Clk de 


flo 


Og flo 




« 2x 10 



-31 



lO'^GeV / V lO'^GeV 



(55) 



The scale factor ratio fl^t/fle is computed directly from the solution of (p with the potential U{(p). The mode function 
Xk is then obtained by solving the equation of motion 

Xk + 3HXk + ^Xk + mlXk = (56) 
fl^ 

with the Bunch-Davies initial condition. 

In Fig.|4] we plot the allowed parameter space given the isocurvature perturbation and relic abundance constraints. 
We see that large local non-Gaussianities by the superheavy dark matter with a small isocurvature power spectrum 
are attainable in the vicinity of the thick dashed line of a = 0.7. Thus, the upper left parameter region of the dashed 
line is ruled out due to overproduction of isocurvature perturbations. Furthermore, the region above of the solid line is 
excluded by the relic abundance condition D.x < iijvi ^ 0.2. These conditions for large local non-Gaussianities yields 
a robust bound on the reheating temperature of 

7kH < 10*^ GeV, (57) 
as well as a bound on the mass of the superheavy dark matter of 

— <2.3, (58) 

which corresponds to mx < A-Hg. We see once again that large observable non-Gaussianities can come from dark matter 
particles that only compose a small fraction of the total dark matter In Fig. |4j this region corresponds to the region 
along the dashed curve that is far below the solid curve. This means that the CMB non-Gaussianities are sensitive 
to bosonic stable particles that are negligible as far as their contribution to the gravitational energy budget today. As 
discussed in Sec. IV the dashed curve in Fig. |4]does not continue indefinitely as nix/m^ since perturbation 



theory for the primordial spectrum breaks down when (% is less than around 10^^. Parametric ally, this breakdown 
point occurs for mx/m^ « 1.4, which corresponds to mx ~ 3/4 . 

From Fig. |4] we can also see that 7= — 1 ~ 0.1 from the isocurvature spectral index nx- This yields the advertised 
result that the isocurvature non-Gaussianities can reach f^i ^ 50. A similar conclusion (in a context different from 
gravitational particle production) was obtained in |69|, where ft^i « 30. 

In this chaotic inflationary model that realizes large non-Gaussianities, the low reheating temperature Trh <C Hg 
implies a long period of matter domination that may lead to nontrivial dark matter and inflaton condensate clustering 
phenomenology (see e.g. rt76l478l ). As demonstrated in Sec. |IV] a lower inflationary scale model that still realizes 
a ^ (9(1) and 7^ 0(0. 1) would allow for an evasion of any phenomenologically undesirable features that may arise 
from the long duration of clustering. 



VI. CONCLUSIONS 



In this paper, we have analyzed the effect of nonthermal dark matter consisting of weakly interacting gravitationally 
produced X bosons on the two and three point functions of the primordial CMB spectrum. We have demonstrated that 
large local non-Gaussianities characterized by an effective /^ri ^ (9(50) can result for a particular set of masses and 
reheating temperatures without violating isocurvature bounds. The conditions that result in large f^i yield a bound on 
the reheating temperature of around 10^ GeV, as well as a bound on the dark matter mass of around 2>Hg < mx ^ AHg. 
For lower allowed values of mx masses, the X bosons can generate a large value despite the fact that they are an 
essentially negligible fraction of the dark matter in the universe. 
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Figure 4: The right bottom of the parameter space under the thick lines is allowed by constraints: the current matter density (the 
solid line) and the upper limit (a ~ 0.7) of the isocurvature power spectrum (the dashed line). The background color shows the 
power spectrum amplitude and the labels on the thin lines represents the spectral index of the isocurvature nx- The right y-axis 
shows the spectral index of the curvature perturbation iij. 



In addition to considering a novel scenario for generating non-Gaussianities, we have have explicitly checked in 
Sec. [B] that the cross-correlation between the entropy perturbation and the curvature perturbation is negligible as is 
required by current data. This computation should be useful for future projects in related topics. Although explicit 
numerical computations were carried out to find the phenomenologically viable region only in the chaotic m^(j)^/2 
inflationary model, we have presented analytic arguments to demonstrate that a similar semi-quantitative behavior 
is expected for most slow-roll inflationary models, including models with lower inflationary scales. The agreement 
between the analytic considerations and the numerical results represents a nontrivial self-consistency check. 

The mechanism presented in this paper connects nonthermal dark matter phenomenology to inflationary phe- 
nomenology. The possibility that a future discovery of large local non-Gaussianities may provide support for the 
existence of a nonthermal dark matter component is indeed intriguing. This may even be considered one of the several 
generic string phenomenological signatures associated with nonthermal dark matter components such as those consid- 
ered in ||79ll80l . although a careful generalization of our work is required to apply the results in that context. We look 
forward to studying further hints nature may be willing to reveal in this direction through observations at the on-going 
Planck mission and other experiments that probe large scale structure. 
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Appendix A: Mode functions beyond the dS approximation 



To calculate the correlation function or the energy density, it is necessary to solve the following mode equation: 

fe2 



Xu{t) + m{t)Xk{t) 



+ mi \Xk{t)=Q. 



(Al) 



In dS space, where H = d/ais constant, this has a well known solution with the Bunch-Davies boundary condition 



Xf(r) 



in which v ~ y'9/4 — m^///^ and //i'' is a Hankel function. In the long-wavelength hmit, k/a <^ H, this solution 
behaves as 



2fl3/2VH 



„;f(v+i/2)^(i) 



k 

oH 



(A2) 



dS 



k 

oH 



r(v) 



(A3) 



in which the decaying solution has been dropped. Fo mx/H > 3/2, \X^^ p o= a^^ dilutes like pressureless dust. 

However, in a quasi-dS spacetime in which H varies slowly, this solution fails to be a good appr oxim ation after 
many efolds past the horizon crossing: At = t — t^, < l/eH, where e = —H/H^. Conversely, Eq. (A2i is a good 
approximation up to the time U just after the horizon crossing. (jATJ can be approximated for f > as 



Xkit) + 3H{t)Xk{t)+mlXkit) « 0, 
and hence we can use the following ansatz for the non-decaying mode for small e <C 1 . 



x,(0=xf(f,) 



a{t) 



exp 



dtiH{ti)v{ti) 



Now consider a quasi-dS spacetime with the constant slow-roll parameter e, e = 0, for which 

1 1 / N 



H{t) Hit,) 
a{t)^a{u){\+eHjfl^. 



Hence, we find 



jjt'H{t')v{t') = i 



v(f)--tanh 



-1 



2v(f) 



- V* 



- tanh 



-1 



3 

2v, 



(A4) 
(A5) 

(A6) 
(A7) 

(A8) 



We note that the mode function Xi^ is oscillatory for imaginary v, but that is not reflected in Eq. ( A8 1 because we have 
kept only the non-decaying mode in Eq. (A3 1. To see the oscillatory behavior of the mode, the decaying mode should 
be taken into account when the real to imaginary transition of v occurs. We thus arrive at 



2-2+2V, 



a{t) 
k 



exp 

-2v. I 



v(f)- - tanh" 



2v{t) 



-V* 



- tanh 



rfv* 



H,a^{t) 



exp 



^5R<( 2v(f)-2v* +31n 



,(v(0 

in which the subscript * indicates the variable is evaluated at and v* is a positive real number. 





^ 3 


'1 




1) 


M*. 




mx 


3^ 


Hit) 


2, 


mx 



(A9) 
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Appendix B: and ^5^^ Correlators 



As shown in Eq. ( 50 1, the non-Gaussianities from this isocurvature model depend on the adiabaticity parameter a 
defined in Eq. ( [48| . Current observational limits constrain the parameter a through the curvature and isocurvature 
correlations |7, 69|. In this section, we will consider the correlations between the curvature perturbations C, and the 
isocurvature perturbations 5s, and show that the correlation is negligible. 

On large scales, the temperature anisotropy is given by 



AT 12 



where in the comoving gauge. 



PCDM PCDM 



= (Ox — — = (Oxds^. 
Px 



(Bl) 



(B2) 



Using the ADM formalism fWT] with the gauge choice 5(p = 0, B = 0, A = 2^, the metric perturbation is written in 
terms of the variable ^. For the metric perturbations defined as Eq. we therefore have 



A 
E 

F 



2C, 
H ' 



aH 



where e = —H/H^. The Einstein equations then yield 



|(ln(«3,))C-5c^O. 



(B3) 
(B4) 

(B5) 
(B6) 



For notational convenience, the superscripts {c){p) will be suppressed in this section. 

To estimate the correlator {t, dx), we have to clarify the possible factors that cause non-zero correlations at leading 
order in (CO- First, since the definition of the energy density given in Eq. ( 18 1 contains the metric g^y, the associated 
perturbations have an explicit ^ dependence, which may introduce non-zero correlations. We can write the energy 
density of X in the comoving gauge as 



Px = 



1 



X' 



2— -rm|x2 



H 



{vxy 



a 



H 



X 



(B7) 



We note that the terms outside the square brackets all decay as a^^, and thus can be neglected at late times. 

The other possibility is due to the effects of metric perturbations on the evolution of X, which appear as interaction 
terms in the Hamiltonian. In terms our gauge fixed variables, the interaction Hamiltonian can be written as 



Hi{t) 




3C U 



mlX^ 



After canonical quantization, the fields ^ and X are mode-decomposed as 



(B8) 




C(r,x) 
Xit,x) 



d^k 



a-jXt{t)e''' + alJCl{t)e 



— ik-x 



(B9) 

(BIO) 
(BID 
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in which the creation and annihilation operators satisfy the usual commutation relations 



q 



(B12) 



The mode functions i^<.(f) satisfies the Bunch-Davies boundary conditions while ^^ (r) satisfies the related adiabatic 
vacuum boundary conditions ll47ll65H67ll . 

In the interaction picture using the in-in formalism, the correlator becomes to leading order 



(B13) 



Here we will restrict ourselves to only few terms in Eq. (B8 1 to estimate (C^x)- Other terms can be estimated in an 
analogous manner Their contribution is similar to or less than what we will estimate here. 

Given that the mode functions Xi^ behave like a^^/^ when tnx > \H and the mass range we are interested in is 
nix ~ He, as well as the fact that its density dilutes as a^^ after the inflation, the density perturbations of X will be 
nearly constant by the end of inflation. Hence we will estimate the density perturbations evaluated at this time. For 
further simplification, we will use px ^ :X^ : /2 and px{t) = {■ Px ■), such that Bs^ ~ : X^ : /2px- 



We first consider the contribution of the interaction term Hj i = / d^xa^3^m^X^ /2 e Hj as follows: 



// dt'{[H,,dt'),C(t,x)5s{t,y)]) =ij dt' j d\{[a\t'%{t\z)-mlx\t\z)X{t,x)(0x5x{tm 



i / dt' / d^xa^{t') 



3 , 3mja)x 



4px(f) 



C(f',z)x2(f',2),C(f,x):X2(fJ): \(B14) 



where we note that the effective coupling constant for the cubic interaction is 2>m\/2. By taking the Fourier transform 
of Eq. (B 14 1, we have 



j d\d'ye-'P-'-'^^-y{i;{tM{t,y))i - {2n)'5^'Hp + q)^A}s,,i{p), 



where subscript 1 denotes the contribution by /// 1. Hence, the power spectrum A^^^ ^ (p) of Eq. (B15 i is 

„3 



Ap) 



3m 



2n^ \Px{t) 



d'k 



(B15) 



(B16) 



With the dS background assumption, the solution of the mode function equation with the Bunch-Davies initial condi- 
tion is 



H 

2Mp^/eVk^ 



(B17) 



2aV^VH 



J{kI2)(v+\I1) 



r(v + i) 



kxY .r(v) f kT 



when-A:T<l, (B18) 



where Hy is a Hankel function, T = f dt\/a{t\) = — \/aH, and v — J 9 / A — m\/Hl? We neglect the contribution of 



t' < in Eq. (B16l because the mode function and Xj^ are oscillatory. Moreover, taking the imaginary part forces 
us to pick at least one decaying mode among the mode functions and Xf; in Eq. (B16 1. Since the decaying mode Q 
and Xi; respectively behave like a^^ and a^^/^^*', the leading imaginary contribution is 



2n^ V Px 



COx 



( dt's [a3(OCp(OCp(0] y ^x,(f')x;(0Xj__~,,(Ox*__jj(0. 



(B19) 



^ Although diverges because of \/e in the exact dS limit, it is unimportant because the computation of interest depends only on the properly 
normalized . 
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The momentum integral has the same form with two-point function of energy density px except the times of two points 
are different. According to Eq. (B18 i, behaves like a^^^^^^, which implies 



2P1 



a{t) 
4?) 



-3+2v 



(B20) 



where is the isocurvature power spectrum frozen after inflation. Setting —kT = 0, we have 



AMjep^ 3 \a{t')H 



(B21) 



where A| — /Sn MpE. Therefore Eq. (B16 1 is simplified as 



ip) 



< 



Px 



a{te) 



-3+2v 



(B22) 



where —3 + 2v < and we have set fg as the time at the end of the inflation because at least by that time, the 
\Xp\^ o= a^^ behavior commences. We note that the \Xp\^ °^ a^^ behavior actually commences earlier than the end of 
inflation because mx/H > 3/2 condition is fulfilled before the end of inflation. However, since we are bounding the 
cross-correlation in this section, we can use this as a conservative estimate. Since px ^ (9(0.01)//'* when mx < H, 
Eq. ( B22 1 can be estimated as 



ip) 



<^Alip)A%ip)('^, 

COx ^ ^ \a{te) 



-3+2v 



(B23) 



As -ln[fl(r,)/fl(f,)] -50 - 70 and A^ - a^'Al 



io-^a2 



< A| if mx/H < 0.7 and cOx ^ 1- 



Eq. ( |B23| l can be used to obtain a clean bound on the cross-correlation and a. We note that 



a{te) 



-3+2v 



\Xp{t,) 



(B24) 



where we assume for the purpose of a conservative estimate that v remains real until the end of inflation as explained 
above. With the assumption that mx — H and Eq. (|34ji, the isocurvature power spectrum in Eqs. p9|) and ([35]) is 



4Ap) 



coi 



\Xp{te) 

\Mt*) 



in which we have dropped a (9(1) factor and we have used 



2a\u.)H{t^)' 



(B25) 



(B26) 



This result intuitively makes sense since the isocurvature fluctuation is — (9(1) at the horizon crossing during inflation 
and then the fluctuation is suppressed by the decreasing mode function Xi^ and the factor cOx until its magnitude 
becomes the order of a'/^A^. Using A|^ (/?) — aA|(p) from Eq. (48 1, we obtain the result 



a{te) 



-3+2v 



(Ox 



aAjip) 



Combining Eqs. (B23 1 and (B27 1, we can conclude that 



(B27) 



(B28) 



unless a ^ A^ 



10 With similar arguments, we can see that the other terms in the interaction Hamiltonian of 
Eq. (|B8|) contribute the same order or less. 
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Appendix C: Justification of tlie background 

In situations in which quantum operators i^^'^' do not have classical expansion field values \ it is necessary 

to justify how the i^q^^ are to be identified. In such cases, the typical procedure is to solve for directly in the 
presence of linear metric fluctuations because quantum operators are generically not spatially homogeneous even in 
Minkowski space while their matrix elements can be. Hence, it is convenient to construct ff^^^ from matrix elements 

of Here we outline how to extract through a spatial average in both the classical and the quantum case. 

In the quantum case, the procedure is to construct a semiclassical quantity that corresponds to the most probable 
semi-classical configuration computed from a quantum expectation value. As explained below, this semiclassical 
construction is meaningful when the quantity has a classical interpretation. 

Let us define the perturbation 5 ff^'''' through ff^''^ = ff^''^ + 5 . If is independent of spatial coordinates for 
some given fixed coordinate choice and 5 ^'''^ satisfies the condition 

lim — / d^x5ff^''^ =0, (CI) 
Vl Jvl 

then we have 

1 f 

df' « Hm — / (f'xe^^y (C2) 

Here we do not use a 3-diffeomorphism invariant measure such that metric perturbations need not be included in the 
spatial average. Since unperturbed quantum operators are intrinsically inhomogeneous, spatial averaging cannot be 
used to extract S''^'' . 

On the other hand, for the quantum case, we are only interested in using it to obtain stochastic boundary conditions. 
For a large class of operators, for any fixed time slice E, there exists a probability functional pzii^e''}) such that 

{d,'"\..e'f'^) = fY\Dd^^pz{{ff'f^})ff^^'\..ff'f'-\ (C3) 

where the left hand side is computed in a fixed vacuum. Since we are computing this at a fixed time, the right hand 

side functional measure Di^'f^ is only discretized over the 3-space of E. It is important to think of an {d^''} element 

in the ensemble (governed by pi) as a classical configuration only for a fixed time slice. The time evolution of {^e'^'} 
may not be governed by classical equations starting from that initial condition. Furthermore, if pz is sharply peaked, 
then the system is essentially in a single configuration (i.e., most elements in the ensemble are similar). This is one 
way of reaching classicality. Suppose there exists a such very probable field configuration on a fixed time slice. It is 



convenient to express this probable configuration (to use for ff'f'' in Eq. (|C2[)) using matrix elements. 



To construct this probable configuration, we start with the following question: if there are samples drawn from 
a quantum governed ensemble, how many would have exactly homogeneous configurations compared to those with 
inhomogeneous configurations. The set of exactly homogeneous configurations form a set of measure zero even though 
it can be the peak of the pz functional (unless this functional diverges at that field configuration point). Most of the 
samples would be inhomogeneous. Having oriented ourselves, the next question that is relevant for us is what is the 



most probable value of Eq. ( C2 1 given samples in the ensemble governed by pz- The number of configurations with 



a given characteristic T is given by 

NY=Nj Y\De'^^^pz{{ff'f^})5{Y)d^i^^, (C4) 

where 5 (F) represents an appropriately generalized delta-function in the functional space and the determinant is there 
for the appropriate normalization. For a fixed spatial average we have the regularized constraint 



* This argument requires generalization when the operators involve time derivatives. This argument will apply for the case of our interest for 
conelators of . 
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in which case A^r gives the number of elements in the sample that realize the homogeneous value of 



^0 



within the volume V^. We will in the end take °°if the infrared regulator removal is meaningful. Hence, we find 

1 



(C6) 



yt(E)' 

which is independent of ff^. The functional derivatives are taken only with respect to functions of 3-spatial variable. 
Incorporating this into Eq. ( C4 1, we find 



(C7) 



Now, the maximum of A^r is obtained by taking a derivative with respect to (f ) (since we are working on a 



single time slice, this derivative need not be functional): 

d 



VVl(E) Mir 



(C8) 
(C9) 



where the prime on the delta-function corresponds to the derivative with respect to the functional argument of the delta 
function. The functional argument of the delta function can be considered to be a function of the variation 



1 



5e\ 



Hence, an integration by parts will yield a solvable equation to the problem of maximizing Nr , . as 



/n 



If ^q"^' is chosen to be ffl'^' such that the ^g'" configurations that satisfy 



1 



i^Ml) Ml) 



d X 



'0 



= 0. 



also satisfy 



we have a solution . Hence, we must look for the peak of 
Consider 



(CIO) 



(Cll) 



(C12) 



(CI 3) 



(C14) 



Assuming pzii^e'^}) is sharply peaked, consider 

/^ln/7z({^i^'}). 

The peak is located at tt^^ = corresponding to the field configuration that satisfies 



(C15) 



(C16) 
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Hence, the quadratic expansion of / about the peak configuration is 

5V 



.f=.f{4'^)+\j 



dx J dx'2 



in) 



4''\^i)) {4"'\x2) - 4'"\^2) 



(C17) 



where the repeated qi indices are summed. ^''^ can be raised in Eq. (|C14| using the usual trick of introducing a source 



5 

^^exp 



7=0 



and carrying out the leading saddle-point approximation integral to obtain 



is) I 



where we used 



(CI 8) 



(C19) 



(C20) 



n 



exp 



0/ 



dx Y dx'2 



5V 



5i^'^''{xi)5ff^^'-\x2) 



We note that if we use a spatially translation invariant state to take the expectation value, ff'^\x) is automatically 
spatially translation invariant and thus 



(C22) 



in Eq. ( |C12| l. Hence, Eqs. ( |C19| l and ( |C22| l combine to give the most probable spatially averaged configuration to be 



.(■0 



(C23) 



to leading order in saddle-point approximation (an expansion in the peakedness of the distribution function pz). 
Hence, when matching to the classical fluid, it is appropriate to consider the homogeneous background quantity 

associated with the quantu m ope rator to be {fff'). It is important to remember that since ffi"'^ (in Eq. (C12 1) is being 

identified with ff^^^ (in Eq. ( C2 1), we are now assuming that the e = solution of perturbation theory is governed by 
equations that depend only on time for the coordinate system that we have chosen. 
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